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ABSTRACT 


Navier-Stokes solutions are used to calculate oscillatory components of pressure, velocity, and 
density, which in turn provide necessary data to compute energy growth factors to determine 
combustion instability. It is shown that wave instabilities are associated with changes in entropy 
and the space and time averages of oscillatory components of pressure, velocity and density, 
together with the mean flow field in the energy equation. Compressible laminar and turbulent 
flows and reacting flows with hydrogen/oxygen combustion are considered in this paper. The 
space shuttle main engine combustion/thrust chamber is used for illustration of the theory. The 
analysis shows that the increase of mean pressure and disturbances consistently results in the 
increase in instability. It is shown that adequate combustion instability analysis requires at least 
third order nonlinearity in energy growth or decay. 

1. INTRODUCTION 

It is well known that growth or decay of energy is responsible for instability or stability of 
waves in fluids, respectively. We recognize that there are three different sources of energy growth 
or decay. They are: (1) pressure fluctuations (acoustic waves); (2) velocity fluctuations 
(hydrodynamic waves due to transition to turbulence, shear layers, or shedding of vortices); and 
(3) density fluctuations (intrinsic waves due to compressibility and/or chain reactions of unstable 
chemical radicals in combustion). Extensive research works have been carried out for acoustic 
instability [1-6] and hydrodynamic instability [6-9], with limited studies available on intrinsic 
instability [10]. Flandro [11] presented the energy balance method in which the acoustic energy 
equation was used to develop the nonlinear stability equation based on isentropic assumption and 
linear superposition of entropy change. Each of the above types of instability requires different 
methods of analyses and unification of such methods has not been attempted. The objective of 
the present paper is directed toward our desire to introduce a unified general method for wave 
instability analyses. 

The basic approach toward this achievement, called the entropy controlled instability (ECI) 
method, is founded on the concept of entropy changes in which nonlinearity is a prime factor. This 
is because our ultimate goal is a general theory which enables nonlinear waves to be studied in 
detail without making simplified initial assumptions such as isentropy. Nonlinear waves invoke 
nonisentropy in which growth or decay of energy is properly accounted for. This line of thinking is 
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dictated also by our desire to deal with both incompressible and compressible fluids, with the 
general approach being capable of handling compressible or reacting flows and a special case 
reduced to incompressible flows. To this end, we examine the energy equation written in terms of 
entropy changes as well as fluctuations of pressure, velocity, and density. Introducing spatial 
averages, or applying the Green— Gauss theorem, we integrate the energy equation by parts 
through which boundary surface integrals arise, playing the role of acoustic intensities. Additional 
terms representing the spatial growth of energy also arise in this process of integration by parts. 
Performing the time averages, we Anally arrive at the nonlinear, nonisen tropic stability equation, 
which assumes the form of nonlinear integro— ordinary differential equation for the energy growth 
factor indicative of stability or instability of waves. The integrands of this nonlinear stability 
equation (the word *nonisentropic > omitted but implied hereafter) are contributed from the results 
of the Navier— Stokes solutions using the Taylor— Galerkin finite elements [12—16]. The nonlinear 
stability equation will then be solved using the fourth order Runge — Kutta method to calculate the 
energy growth factor. It is shown that the general theory presented here is reduced to a simple 
linear theory of Cantrell and Hart [1] or Culick [2] as a special case with appropriate assumptions. 

We shall discuss the details of these processes in the following sections. 

2. ENTROPY CONTROLLED INSTABILITY (ECI) METHOD 


2.1 Navier — Stokes Solutions 


In order to perform stability analysis, it is necessary to first solve the unsteady Navier— Stokes 
equations representing the desired physics. To this end, the most general conservation form of 
Navier Stokes equations for compressible flows is written as 
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where Ty is the viscous stress tensor 
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and E is the stagnation energy 



E = e + i v iVi = c p T-H + ivi v i (If) 

and is the body force and qj is the heat flux vector. 

N 

qj = - AT,j + P D ^ H k Y k ,j (lg) 

Here, A and D are the thermal conductivity and mass diffusivity, respectively. H k is the total 
enthalpy of species k, Y k is the mass fraction for the species k, and w k is the reaction rate for the 
species k. Note that all derivatives are covariant in case of cylindrical coordinates. 

For turbulent flows we add to (1) additional transport equations for turbulent kinetic energy 
and dissipation energy (k — € model), respectively, 

IM + hfn k -tw) = G ~ pe ( lh ) 

+ ~ = c i i G ~ c 2^ f2 /k (ii) 

where G is the turbulent thermal dissipation energy transport and 
**t H k 2 

* = #I+ V * = * + * as ' c n7 * 

Cji = 0.09, _ c^ = 1.44, c 2 = 1.92, cr k = 1.0, <r € = 1.3 

Appropriate modifications to the equations of momentum and energy and the suitable turbulent 
combustion model should be included as detailed in [12]. 

The solution procedure using the Taylor—Galerkin method are found in [13—17], which will 
not be repeated here. The space shuttle main engine combust ion/ thrust chamber is used for 
numerical analysis. Examples include shock waves interacting with laminar and turbulent flows 
and with laminar reacting flows. Calculations of energy growth factors and stability criteria are 
discussed below. 


2.2 Nonlinear. Nonisentropic Stability Analysis 


With the solutions of unsteady Navier— Stokes equations for the time— dependent periodic 
oscillatory initial boundary conditions available, we now turn to the entropy controlled instability 
(ECI) method of determining the unstable wave phenomena. To introduce entropy changes in the 
energy equation, we begin with thermodynamic relations for an ideal nonisentropic gas and write 
the pressure gradients in terms of density and entropy gradients, 

2 

P>i= a *i + — S -i (2) 

C P 

where a is the speed of sound, S is the specific entropy, and the comma denotes the partial 
derivative with respect to x^ The gradient of the stagnation energy, E, assumes the form 


E, i = ( C P T- f + * v j v j)'i 


where the repeated indices imply summing. Performing the differentiation implied in (2) results in 
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( 4 ) 


E -i = R^ P ’i-R- V*i + v j v i‘i 

P 

where R is the specific gas constant. Substituting (2) into (4) yields 

pE,i = H p ti + g S ti + p v jVj)i 

Let us now examine the conservation form of the energy equation 

Jt M + ( pE v i “ a U v j)’i = 0 

where a y is the stress tensor 

a ij = p*ij + ( v i>j + v j-i -| v k>k tf ij) 

Substituting (5) into (6) gives 

It W + E ( p v i)*i + v i[? ^’i+R S ’i +pv j v j>i]-( a ij v j)-i = 0 

This is the entropy controlled energy equation, instrumental in determining the nonlinear 
instability. 


The pressure, velocity, and density may be written in the form 
P = P + p' 

V 1 = + y[ 

p = p + p 


where the symbols, bar and prime, denote the mean and fluctuation parts, respectively, 
thermodynamic relations we may write the entropy difference in the form 


1 7 



Expanding the RHS of (10) in infinite series and after some algebra, we obtain 
S = R [S( t) + S ( 2) + S ( 3 ) + S ( 4) + ...J + S„ 

where S 0 is the entropy at the initial state (S Q = 0), and 
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Here the terms with fifth order or higher are neglected, assuming that they are negligible. 

To obtain acoustic energy equation, we apply the Green— Gauss theorem or integrate (8) by 
parts and take a time average in the form 

<J n ftW dft + J r [ E ^ v i n i + v i(p n i + ! s n i + p v j v j n i) n i] dr - J n [ E *i t"\ 

+ P ( v i f)'i + ( v i f)*i S + (^iVj).i vjdft ) = 0 (12) 

where fi and T represent the domain and boundary surface, respectively, ( - ) implies time 
averages and n^ denotes the component of a vector normal to the surface. It is interesting to note 

that the domain integral terms (last four terms) arise as a result of the integration by parts 
implying the spatial growth of energy. This spatial growth of energy is in addition to the usual 
temporal growth of energy indicated by the domain integral with the time derivative (first term in 
(8)). Let e be the energy growth factor, stable for 0 < € < 1, unstable for e > 1, with c=l 
indicating the neutral stability. The fluctuation terms in (9) and (11) contribute to (12) as 
multiples of higher order fluctuations. Multiplying the like powers of € with the fluctuation terms 
of the same order and retaining the terms up to and including the fourth order, we obtain 

<J(0 dn> = <J (* 0 + ef, + e\ + + e 4 * 4 + ...) dft) (13a) 

(J (*) dr) = + f<t> , + « 2 0 2 + f3 ^3 + + — ) dr) (13b) 

Here S Q and contain only the mean quantity, S 2 , £ 3 , S 4 and <j > 2 , ^ 3> ^4 are referred to as 

the first, second, third, and fourth order perturbations capable of growth and decay in acoustic 
energy as dictated by the magnitudes of e, growing if e > 1 and decaying if e < 1. Notice that 
the relations in (9a, b, and c) imply that the Navier— Stokes solutions are obtained with the initial 
condition 6 = 1 and the energy in (13) may grow or decay in accordance with the magnitude of e 
from the initial or reference state e = 1. 


Performing the differentiation as implied in (12), we obtain 

If G 2 Et + e 3 E 2 + e 4 E 3 ) = t \ + e 3 I 2 + e 4 I 3 (14) 

in which the zeroth order terms are canceled and first order terms vanish due to time averages. 
Notice that the time averages are contained in E( ^ and I( ^ . Here, the energy growth factor e is 

an explicit function of time. However, E( £j is no longer an explicit function of time because of its 

time averages. Thus the partial derivative with respect to time in (14) involves only e, not E( ^ , 

so that 
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or 


3E. 


where higher order terms and those terms much smaller than unity have been neglected. 
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It follows from (15) that the nonlinear stability equation takes the form 

- a,* - a 2 « 2 - a 3 «* = 0 (16) 

This is a form identical to Flandro [11], although the basic approach to the formulation and the 
solution procedures differ. Here a lf a 2 , and a 3 are energy growth rate parameters of first, second, 


and third order, respectively. 
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Explicit forms of ^ , and C ( ^ will be discussed 

in Section 3 and the integrands for the 


second and third order growth parameters are shown in Appendix A. The physical significance 
and interpretations of the above process will be discussed in the following subsection. 


2.3 Physical Interpretation 

The physical significance in the above development may be schematically demonstrated as 
shown in Fig. 1. First, unsteady Navier—Stokes solutions with time dependent oscillatory initial 
boundary conditions are obtained at each computational grid as shown in Fig. la. A typical 
variable vs time, (Fig. lb), at any grid point (finite element node) may exhibit sawtooth type 
wave forms as well as sinusoidal waves. Imagine that several variables at each of the several 
thousand nodes present themselves in the form as shown in Fig. lb. Each wave form of different 
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frequencies at different spatial locations is expected to contribute to the overall stability or 
instability. Take one complete peak wave form period as indicated by r = nAt where At is the 
Navier— Stokes time step with n usually being in the range between teens and hundreds. 

Performing the time average for this time period (nAt), we obtain the mean quantity (f) of 
the variable f. Then the disturbance (fluctuation) part f ' is calculated as 

f ' = f- 7 

where f is the Navier— Stokes solution. 

From these data, spatial integrations involving both domain and boundary surfaces and time 
averages (Fig. la) as dictated by (18) and (19) are then carried out. 

Here, a (n) dfl denotes the temporal growth of energy, whereas J 6 (n) dfl represents the 

spatial growth of energy. The boundary integral J c[ n) n £ dr indicates the acoustic intensity. 

For combustion chamber applications, the burning surface admittance, response functions, or 
information on nozzle entrance or acoustic liner, etc., can be initially imposed on the 
Navier— Stokes solutions as boundary conditions to generate the oscillatory motions as depicted in 

Fig. lb. These boundary data, however, reappear as called for in the integrand, C^ n) , in (19). 

They act as acoustic intensity and eventually contribute to the determination of stability or 
instability of the system. 


The ingredients of a (n) , 6 (n) , and c[ n) arise from fluctuation parts of the Navier— Stokes 

solutions as well as the mean parts. These fluctuations through the temporal growth of energy 

a (n) , spatial growth of energy b ( n) , and acoustic intensity c[ n) , are responsible for driving the 

system toward instability. The consequences lead to determination of the energy growth factor vs 
time as shown in Fig. lc by solving the nonlinear stability equation (16). It is important to realize 
that the Navier — Stokes solutions and the nonlinear stability equation (16) encompass acoustic 
waves, hydrodynamic (vortical or shear layer) waves, intrinsic (density fluctuations or chemically 
reacting combustion) waves, or combined effects of all types of wave interactions. 

Notice that, in (6) and (12), heat flux changes are not directly involved. This is because the 
spatial integration of fluctuations of heat flux or temperature vanishes upon time averages. 
However, fluctuations of temperature and chemical species mass fractions have been included in 
the Navier— Stokes solutions of (1), thus indirectly contributing to the fluctuations of density (p) 
in the stability calculations. 

Some comments on the mean pressure changes (DC shifts), pressure coupling, velocity 
coupling, limit cycles, and triggering are in order. These phenomena often observed in solid 
propellant combustion chambers, would appear in the flowfield under appropriate initial and 
boundary conditions during the Navier— Stokes calculations. The purpose of the ECI method is 
merely to determine whether the flowfield containing such physical phenomena is stable or 
unstable. It is emphasized that all physical phenomena are expected to prevail and be identified 
in the solution of Navier— Stokes equations, which are then reconfirmed in the stability analysis. 


3. LINEAR INSTABILITY 

To gain an insight into the solution of (16), we may neglect the last two terms associated 
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(20) 


with the second and third order energy growth rate parameters and write 

which yields a solution in the form 

in € = fl t t + c t (21) 

To establish an initial condition, we assume neutral stability € = 1 at t = 0. This gives c 1 = 0. 
Thus, the solution takes the form 
a.t 

* = e 1 (22) 

Under this initial condition, there exists a unique solution for any given a ^ with t > 0. The 
stability criteria are: 

Stable: 0 < e < 1 with — ao < < 0 

Neutral stability: e = 1 with = 0 
Unstable: 1 < € < m with 0 < < m 

Although these criteria are not applicable for the nonlinear equation (16) when fl 2 anc * are 
involved, the same initial condition, = 1 at t = 0, as originally defined in (13a, b), can be used. 
That is, there exists a unique solution e for any given fl lf a 2 and a 3 with t > 0. Therefore, the 
criteria for stability in terms of e remain the same with various combinations of a v fl 2 , and fl 3 . 


3.1 Nonisentropic case 


Notice that the integration involved in the growth rate parameters fl^, a 2 , and fl 3 and the 

subsequent solution of the nonlinear ordinary differential equation (16) are formidable. However, 
it would be informative to examine the case of linear instability (a 2 = 0, fl 3 = 0) given by (20) 

but with nonisentropy. Therefore, the first order energy growth rate parameter takes the form, 
from (20) and (17a), 
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Assuming that the fluid is inviscid (/* = 0) the explicit forms of integrals in (23) become 
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Notice that the first term on RHS of (24a) represents the kinetic energy of the sound wave. 

If isentropic assumption is made, then p can be expanded into infinite series in terms of p’, p’ 2 , ... 
such that the second term on the RHS of (24a) contains the potential energy of the sound wave 
usually identified in the linear acoustic energy equation. However, our objective here is to allow 
entropy to change and, therefore, we must keep p to remain subjected to nonisentropy. Another 

important observation is the domain integral for b ( which has appeared for the first time as a 
result of the integration by parts (24b), signifying the spatial growth of energy. The boundary 

integral for c^nj contains the terms of acoustic intensity normally identified in the linear 

acoustics, but significantly in a different form. Additional terms which arise in the process 
developed in this analysis will allow determination of stability or instability of wave motions with 
explicit changes of entropy contributing to the growth or decay of energy. However, it is not 
possible to evaluate the integrals (24a, b, c) because the density fluctuations p cannot analytically 
be determined. This difficulty can be resolved in a certain special case if isentropy is assumed as 
discussed below. 


3.2 Isentropic case 


If the flow is isentropic then it can be shown that 
-P^iplrlP’ 2 | 1 ?(1— 7)(1~ 27 ) P’ 3 , 1 ?(1— 7)(1— 27)(1— 37) p' 4 , 
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where a denotes the speed of sound without flow. As far as the linear stability is concerned only 
the first two terms on the RHS of (25) will contribute to (24a, b, c) as seen from substitution of 
(25) into (9c) and subsequently to (24a, b,c). 
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Since p does not appear in the integrals of (26a, b, c) it is now possible to substitute 
analytical forms of p' and v* in terms of time dependent acoustic eigenfunctions and perform 


explicit analytical integrations. However, it is clear that the additional domain integral (26b) and 
many more additional terms arising as a result of the entropy controlled energy equation given by 
(8) or (12) will produce the results quite contrary to the traditional isentropic solution, (see 
Cantrell and Hart [1] or Culick [2]) 



(27) 


Note that the terms on the numerator are the same as the first four terms of (26c) with the 
negative sign as indicated in (19a), and the terms of the denominator are identical to those in 

(26a). However, the domain integral of (26b) is absent. This integral contributed by b* ^ 
represents the spatial growth or decay of energy as balanced by the boundary conditions and 
acoustic intensities on the boundary surface. This is in contrast to the temporal growth or decay 
of energy originating from (26a). 


Despite the special feature in the proposed formulation, however, the analytical forms for p' 
and v^ as used by Cantrell and Hart [1], are incapable of simulating sawtooth type shock waves. 

For example consider the acoustic field of a cylinder with the radius R and length L, 


p' = 9*e |p cos ^ J m g) cos (m0) 

V ; = 9k [ p cos ^ J m (^ cos (mfl) eiw»l 

L f.tn J 


(28a) 

(28b) 


up 

Here angular frequency, u, is defined as 
»*[(£)’♦(£)]* 


with m, n, being the positive integers or zero, J mn the Bessel function of order m, and ^ the 

nth root of Jm n (/3) = 0. Retaining only those terms arising from the standard acoustic energy 

equations (without employing the entropy controlled energy equation), a simple solution can be 
obtained using (28a, b) to calculate the linear energy growth rate parameter as demonstrated 

by Cantrell and Hart [1]. If the entire terms implied in (26a, b, c) are used, however, it is no 
longer possible to obtain analytical solutions. 
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4. NONLINEAR, NONISENTROPIC INSTABILITY 


Nonlinear, nonisentropic waves occur in many industrial propulsion combustion systems. 
Shock waves may interact with turbulent vortical waves or shear layers of liquid jets in gas 
medium. Density fluctuations due to chemical reactions may also be combined. We have seen 
that analytical solutions of even the linear instability by ECI method presented in the previous 
section are intractable. First of all, the fluctuation variables p 1 , v^, and p are to be numerically 

calculated by solving the Navier— Stokes equations. Then we must perform numerical integrations 
required to evaluate the energy growth rate parameters a lf a 2 , and 0 3 - Subsequently, numerical 

solutions of the nonlinear ordinary differential equation must be carried out. 

Numerical integrations as required by (18a, b, c) and (19a, b, c) can be performed most 
efficiently by Galerkin finite element techniques to calculate the energy growth rate parameters a^, 

a 2 , and a 3 according to (17a, b, c). Finally the nonlinear ordinally differential equation (16) is 

solved using the Newton— Raphson method to determine the energy growth factor e. Thus, it is 
seen that the determination of stability (0 < e < 1), instability (c > 1), or neutral stability (e = 
1) is made available for each nAt period and we move on until desired time is reached as shown in 
Fig. lc. It is interesting to see that, in this decision making process of stability or instability, "all 11 
nodal points (thousands of nodes) with their nodal values of "all" variables have participated. 
Every wave peak, whether sinusoidal or sawtooth type, has been recognized. Shock waves 
interacting with turbulence, shear layers, or shedding of vortices, or effects of chemical reactions 
can be reflected in calculations of (17 a, b, c) and could have eventually contributed to the 
solution of (16) for determination of stability or instability. 

Solutions to the nonlinear ordinary differential equation (16) may be obtained using the 
fourth— order Runge— Kutta method. Iterations will continue until convergence. 

5. SOLUTION PROCEDURE 

It is clear that the solution consists of three parts. First the Navier— Stokes equations are 
solved. Then the results of the Navier— Stokes solutions are used to calculate the energy growth 
rate parameters. Finally the energy growth factor is computed by solving the nonlinear, 
nonisentropic stability equation. Step— by— step solution procedures are described as follows: 

(1) With appropriate boundary and initial conditions, solve the Navier— Stokes equations. 

Initially, the mean pressure, p, and temperature, T, based on the ideal gas law, are specified 

everywhere. At the inlet, however, the oscillatory pressure is input \p = p (1 + d sin art)], 
where d is the % disturbance and w is the frequency. 

(2) Calculate p, v if />, and T. Taylor— Galerkin finite element method with adaptive meshes is 
used in Navier— Stokes solutions. 

(3) Advance time steps (At) of Navier— Stokes solutions to obtain wave oscillations to cover at 

least one wave period. At is determined continuously which satisfies an acceptable Courant 

number. 

(4) Take time averages for the period nAt with n chosen such that at least one peak wave 
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period is covered. These time averages lead to p, vj, and p. 

(5) Calculate the fluctuation quantities as p* = p — p, v| = vj — vj, p s p — p, where p, v^, 
and p represent Navier— Stokes solutions. 

(6) Calculate the energy growth rate parameters atj, a 2 , and a 3 from (17a, b, c) using the 
results of step 5, above. 

(7) Solve the nonlinear differential equation (16) using the Newton— Raphson method with the 
initial condition e = 1 at t = 0, corresponding to neutral stability. 

(8) Repeat steps 1 through 7 until the desired length of time has been advanced. 

Note that for each time— average period in step 4, above, instability and stability are 
determined by e > 1 and < < 1, respectively, with e = 1 being the neutral stability. If the system 
is found to be unstable, it is not necessary to proceed to the next time step. However, for the 
entire ranges of time for which Navier— Stokes solutions are available, the stability analysis may be 
performed, even if instability has been found in previous time steps. This is so because 
Navier— Stokes solutions are independent of the stability analysis as formulated here. Rather, the 
stability analysis here determines the state of stability or instability based on current flowfield as 
calculated from the Navier — Stokes solution. The following applications are based on computer 
code ECI— 2. 


6. APPLICATIONS 


To illustrate the theory, an axisymmetric geometry of combustion/thrust chamber is 
investigated for various types of flow: Case 1, laminar compressible nonreacting flows; Case 2, 
turbulent compressible nonreacting flows, and Case 3, chemically reacting laminar compressible 
flows with hydrogen/oxygen combustion. 

Case 1 L aminar Compressible Nonreactinv Plows 

Figure 2a shows the axisymmetric geometry of combustion/thrust chamber with 965 
triangular elements and 528 nodes as a consequence of adaptive mesh process. Notice that in the 
vicinity of the walls approximately 80% of the nodes are concentrated resulting in prominent 
boundary layers. To demonstrate the capability of the code, initially, the steady state flow field 
without disturbances is examined. The velocity vectors and contours of Mach number, pressure, 
and temperature are shown in Fig. 2b, c, d, and e, respectively. Formation of boundary layers 
(Fig. 2b), separation of boundary layers and weak shock waves downstream (Fig. 2c, d), and 
decrease of temperature downstream and toward the wall (Fig. 2e) are evident. 

With disturbances of d = 10%, 20% and 30% imposed on the mean pressures at the inlet, the 
Navier— Stokes transient analyses are performed. The time steps are continuously adjusted to 
satisfy acceptable Courant numbers, 0.2 < CN < 0.4, for convergence. The graphical 
representation of oscillations of all variables at every node versus time is overwhelmingly complex. 

Therefore, wave forms only for d = 30% and p = 3000 psi at three selected positions, A at (1.8, 
11.43 cm), B at (31.75, 11.43 cm), and C at (63.5, 6.55 cm) are shown in Fig. 3. Note that, 
although the sinusoidal input is provided at the inlet, the oscillations downstream become 
nonlinear, possibly of sawtooth type. 



In Fig. 4, the energy growth factors e for p = 500 psi are shown for various % disturbances. 
It is seen that for d = 10%, the energy growth factor remains in the stable region e < 1. The 
dotted lines and solid lines indicate the results for the linear (Eq. 20) and nonlinear (16) cases, 
respectively. As the disturbance increases (d = 20%) the energy growth factor reaches the neutral 
stability, e = 1 at t 2 0.015 sec. For d = 30%, the energy growth factor increases further (e = 
1.038) at t £ 0.2 sec. Notice that the linear analysis underestimates the stability if stable whereas 
it underestimates the instability if unstable. The general trend is that instability is proportional to 
the percent disturbances. 


As the mean pressure increases (p = 3000 psi), the possibility of instability increases as shown 
in Fig. 5, with the conclusion that instability is proportional to the mean pressure. Recall that for 

d = 10%, p = 500 psi, stability prevailed throughout whereas with d = 10%, p = 3000 psi neutral 

stability has been reached. The peak values of e for p = 3000 psi are significantly larger than 

those for p = 500 psi. 

Case 2 Turbulent Compressible Nonreacting Flows 

The discretized geometry for a steady state turbulent compressible nonreacting flow is shown 
in Fig. 6a with a total of 2688 elements and 1416 nodes. As seen in Fig. 6b (velocity vectors) and 
Fig. 6c (Mach number contours), the boundary layers are thinner than in the case of laminar flow, 
leading to turbulent shock wave interactions. It is clear that gradients of pressure (Fig. 6d) and 
temperature (Fig. 6e) are larger than in the case of laminar flow. 


Wave forms of transient turbulent nonreacting flow with disturbances, d = 30%, p = 3000 
psi, at the three positions are shown in Fig. 7. Although the wave forms in turbulence at these 
positions do not seem much different from the case of laminar flow, it is quite possible that wave 
forms at other locations where turbulent velocity gradients are significant would be drastically 
changed in contrast to the laminar flow. 


The energy growth factors for turbulent flow are shown in Fig. 8 for p = 500 psi and Fig. 9 

for p = 3000 psi. It is interesting to see that the effect of turbulence is to increase instability as 
compared with the laminar flow. The general trend other than the turbulent effect, however, 
remains the same as the laminar flow. That is, instability is proportional to disturbances and the 
mean pressure. The linear analysis again shows that stability is underestimated if stable and 
instability is underestimated if unstable. 

Case 3 Chemically Reacting Laminar Compressible Flows with Hydrogen/Oxygen Combustion 

Figure 10 shows the discretized geometry for a steady state chemically reacting flow without 
disturbances. A total of 1580 elements and 844 nodes are used. The chemical reactions considered 
are shown in Appendix B. The mass fractions at the inlet are 0.111 for H 2 and 0.889 for 0 2 , and 

the inlet velocity is 500 m/s. In this analysis the effect of viscosity is ignored to ensure an 
enhanced computational convergence, which leads to a flow without boundary layers along the 
wall. Due to the finite rate chemistry and stiffness arising from the chemical source terms, the 
computational convergence is rather slow. In the region of chemical reactions the velocity is 
decreased (Fig. 10b) and the Mach number contours are spaced widely apart (Fig. 10c), resulting 
in a rapid increase of pressure (Fig. lOd), but no evidence of shock discontinuities (Fig. 10c, d) 
between the throat and the downstream nozzle area. Temperature increases in the region of 
combustion in a sharp contrast to the nonreacting cases (Fig. 2e and Fig. 6e). 
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The contours of 8 species are shown in Fig. 11, beginning with the reactants, hydrogen (Fig. 
11a) and oxygen (Fig. lib), followed by products, H (Fig. 11c), H0 2 (Fig. lid), H 2 0 (Fig. lie), 

H 2 0 2 (Fig. Ilf), O (Fig. llg), and OH (Fig. llh). It is observed that H 2 is depleted halfway 

between the inlet and throat whereas 0 2 prevails further downstream before it is depleted behind 

the throat. Most of the products (including the radicals), H, H0 2 , 0, and OH rapidly increase 

from the inlet and become maximum as they pass through the throat. In contrast, H 2 0 and H 2 0 2 

gradually increase downstream with a constant rate. 

Figure 12 shows the wave forms of the transient chemically reacting flow with disturbances, d 

= 30% and p = 500 psi, at the three positions considered earlier. Notice that, with chemical 
reactions, the frequencies of waves are very low and the response is quite slow particularly at 
downstream locations. Once again, the results plotted in Fig. 12 are misleading because 
oscillations in other locations (over 800 nodes) may prove to be significantly different. After all, 
combustion instability is determined by the stability equation, not by the appearance of 
oscillations observed at random locations. 

In Fig. 13, the energy growth factors for the chemically reacting flow are examined. First of 
all, the linear analysis (dotted line) shows an apparent faulty prediction of instability for d = 10%, 
in which the nonlinear analysis indicates 6 £ 0 throughout the time segment investigated. 
Obviously, this is an indication of unreliability of the linear analysis. For d = 20%, however, the 
linear analysis appears to give the consistent results similar to the earlier examples in that the 
linear analysis underestimates the stability when stable and underestimates the instability when 
unstable. The nonlinear analysis shows that, for d = 20% and d = 30%, the peak values of energy 
growth factors are significantly larger than the nonreacting cases. Does this imply that chemically 
reacting flows always tend toward instability? An affirmative answer to this question is 
premature. In fact the entire investigation presented here is subject to the future verification by 
experimental measurements. 


7. CONCLUSIONS 

The entropy controlled instability method has been applied to various problems in laminar 
flows, turbulent flows, and reacting flows for determination of stability conditions. The following 
conclusions are reached: 

(1) Instability increases with an increase of disturbances. 

(2) Instability increases with an increase of the mean pressure. 

(3) Instability increases as laminar flows are changed to turbulent flows. 

(4) Instability due to production of radicals under the finite rate chemistry is significant for 
the case investigated. 

(5) The linear stability analysis underestimates stability if stable and underestimates 
instability if unstable. 

(6) The correct stability analysis calls for at least the third order nonlinearity. 

(7) Future studies are required to validate the present theory in comparison with 
experimental results for those cases other than examined in this paper. 
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APPENDIX B 


REACTION EQUATIONS 


H 2 + 0 2 = 20H 
H + 0 2 = OH + 0 
OH + H 2 — H 2 0 + H 
O + H 2 = OH + H 
20H = H 2 0 + O 
H + OH + M = H 2 0 + M 
2H + M = H 2 + M 
H + 0 2 + M = H0 2 + M 


H0 2 + H = H 2 + 0 2 
H0 2 + H = 20H 
H0 2 + O = OH + o 2 
2H0 2 = H 2 0 2 + 0 2 
H0 2 + H 2 “ H 2 0 2 + H 
H 2 0 2 + OH = H0 2 + H 2 0 
H 2 0 2 + H =: OH + H0 2 
H 2 0 2 + O = OH + H0 2 


H0 2 + OH = H 2 0 + 0 2 


H 2 0 2 + M = 20H + M 


Any variable (f) 




(b) Navier— Stokes solutions with initial oscillatory boundary conditions 



(c) Energy growth factor ( e ) 


Pig. 1 Schematic overview of ECI 
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Fig. 2 


r 



88.9 cm 

(a) Discretised geometry, 965 element*, 528 nodes 



(b) Velocity field, 1 mm = 105 m/s 



(c) Mach number contoun, max = 1.6, min = 1.6, incr = 0.032 



(d) Pressure contours, max = 3000 psi, min = 230 psi, incr = 53 psi 



Steady state laminar nonieactiiig flow without disturbances. M (inlet) = 0.2, Re 

= 66,000, /* = 0.23 kg/m 2 /s, 7 = 1-2, p = 3000 psi, p = 55.8 kg/m 1 , T = 
3656.33' K 
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Fir 3 Waveforms for transient laminar nonreacting flow with disturbances. Pic— aic, 

6 density, axial velocity, and radial velocity vs time at various positions. Position 
A (1.8, 11.43 cm), position B (31.75, 11.45 cm), position C (63.5, 6.55 an) 
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Fig. 4 


Energy growth factor vs time, transient laminar flow with disturbances, p = 500 
pgi, f (inlet) = 1022 Hz, M (inlet) = 0.2 
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Fig. 5 Energy growth factor vs time, transient laminar flow with disturbances o = 3000 
psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 ’ P 
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(a) Diacre tised geometry, 2688 elements, 1416 nodes 




(c) Mach number contours, max = 2.2, min = 0, incr = 0.044 




(e) Temperature contours, max = 3656.33* K, min = 15 . 45 * K, incr = 72.82* K 



Fig. 6 Steady state turbulent nonreacting flow without disturbances. M (inlet) = 0.2, 

Re = 500,000, /i = 1.737 kg/m 2 /s, 7 = 1.2, p = 3000 psi, p = 55.8 kg/m 3 , T = 
3656.33* K 
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Fig. 7 


fo * t /i ran 5. ei ; t ^rbnlent flow with disturbances, pressure, density, 

U43 m^Porit^B tU Tv V ^ ri0UiJ P°®‘ioiis. Position A (L8, 

11.43 cm}, Position B (31.75, 11.45 cm), Position C (63.5, 6.55 cm) 
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Fig. 8 Energy growth factor ys time, transient turbulent flow with disturbances, 
p = 500 psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 


0.90606 1.09024 0.92401 1.04672 0.09747 1.00044 





Pig. 9 Energy growth factor vs time, transient turbulent flow with disturbances, 
p = 3000 psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 
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(a) Diacre tised geometry, 1580 elements, 844 nodes 




(d) Pressure contours, max = 1897 psi, min = 46 pel, incr = 37 psi 



Fig 10 


Steady state chemically reacting flow without disturbances, inlet velocity = 500 

m/s, Re = 66,000, p = 0, 7 = 1.413, p = 500 psi, p = 40.08 kg/m 1 , T = 500* K, 
inlet mass fraction H 5 = 0.111, O, = 0.889 
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(a) H, contour*, mu = 0.111, min = 0.11, incr = 0.8 * 10* 5 




(c) H contour*, max = 0.145 * 10' 14 , min = 0, incr = 2.89 * 10' 1# 



Fig 11 Distributions of various chemical species at steady state without disturbances 
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(0 HjOj contoun, max = 0.412 « 10* 3 , min = 0, incr = 0.82 » 10 s 




Fig 11 Distributions of various chemical species at steady state without disturbances, 
continued. 
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position C 


Fig 12 


Wire forms for transient chemically reacting flow with disturbances i 
density, axial velocity and radial velocity vs time at various positions m 
(1.8, 11.43 cm), Position B (31.75, 11.45 cm), Position C (63.5, 6.55 cm) 
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Fig 13 Energy growth factor vs time, transient chemically reacting laminar flow with 
disturbances, p = 500 psi, f (inlet) = 5583 Hz, inlet velocity = 500 m/s 


